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Abstract 

The structure and dynamics of the solar corona is dominated by the magnetic field. In most 
areas in the corona magnetic forces are so dominant that all non-magnetic forces like plasma 
pressure gradient and gravity can be neglected in the lowest order. This model assumption is 
called the force-free field assumption, as the Lorentz force vanishes. This can be obtained by 
either vanishing electric currents (leading to potential fields) or the currents are co-aligned with 
the magnetic field lines. First we discuss a mathematically simpler approach that the magnetic 
field and currents are proportional with one global constant, the so-called linear force-free field 
approximation. In the generic case, however, the relation between magnetic fields and electric 
currents is nonlinear and analytic solutions have been only found for special cases, like ID 
or 2D configurations. For constructing realistic nonlinear force-free coronal magnetic field 
models in 3D, sophisticated numerical computations are required and boundary conditions 
must be obtained from measurements of the magnetic field vector in the solar photosphere. 
This approach is currently of large interests, as accurate measurements of the photospheric 
field become available from ground-based (for example SOLIS) and space-born (for example 
Hinode and SDO) instruments. If we can obtain accurate force-free coronal magnetic field 
models we can calculate the free magnetic energy in the corona, a quantity which is important 
for the prediction of flares and coronal mass ejections. Knowledge of the 3D structure of 
magnetic field lines also help us to interpret other coronal observations, e.g., EUV-images of 
the radiating coronal plasma. 



1 Introduction 



The magnetic activity of the Sun has a high impact on Earth. As illustrated in Figure 1, large 
coronal eruptions like flares and coronal mass ejections can influence the Earth's magnetosphere 
where they trigger magnetic storms and cause aurorae. These coronal eruptions have also harm- 
ful effects like disturbances in communication systems, damages on satellites, power cutoffs, and 
unshielded astronauts are in danger of life-threatening radiation. 1 The origin of these eruptive 
phenomena in the solar corona is related to the coronal magnetic field as magnetic forces dominate 
over other forces (like pressure gradient and gravity) in the corona. The magnetic field, created by 
the solar dynamo, couples the solar interior with the Sun's surface and atmosphere. Reliable high 
accuracy magnetic field measurements are only available in the photosphere. These measurements, 
called vector magnetograms, provide the magnetic field vector in the photosphere. 




Figure 1: Magnetic forces play a key role in solar storms that can impact 
Earth's magnetic shield (magnetosphere) and create colorful aurora. [Source: 
http://sohowww.nascom.nasa.gov/gallery/images/magnetic_clean.html. Courtesy of SOHO 
consortium. SOHO is a project of international cooperation between ESA and NASA.] 

To get insights regarding the structure of the coronal magnetic field we have to compute 3D mag- 
netic field models, which use the measured photospheric magnetic field as the boundary condition. 

1 For an animation of a coronal mass ejection (CME) causing a substorm and aurora, see 
http:/ /sohowww. nascom.nasa.gov/gallery/Movies/recon/reconsm.mpg 
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This procedure is often called "extrapolation of the coronal magnetic field from the photosphere." 
In the solar corona the thermal conductivity is much higher parallel than perpendicular to the field 
so that field lines may become visible by the emission at appropriate temperatures. This makes 
in some sense magnetic field lines visible and allows us to test coronal magnetic field models. In 
such tests 2D projection of the computed 3D magnetic field lines are compared with plasma loops 
seen in coronal images. This mainly qualitative comparison cannot guarantee that the computed 
coronal magnetic field model and derived quantities, like the magnetic energy, are accurate. Coro- 
nal magnetic field lines which are in reasonable agreement with coronal images are, however, more 
likely to reproduce the true nature of the coronal magnetic field. 




Figure 2: Plasma (3 model over active regions. The shaded area corresponds to magnetic fields 
originating from a sunspot region with 2500 G and a plage region with 150 G. The left and right 
boundaries of the shaded area are related to umbra and plage magnetic field models, respectively. 
Atmospheric regions magnetically connected to high magnetic field strength areas in the photo- 
sphere naturally have a lower plasma (3. The original figure was published as Figure 3 in Gary 
(2001). 

To model the coronal magnetic field B we have to introduce some assumptions. It is therefore 
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necessary to get some a priori insights regarding the physics of the solar corona. An important 
quantity is the plasma (3 value, a dimensionless number which is defined as the ratio between the 
plasma pressure p and the magnetic pressure, 

^ = 2mo^- (1) 

Figure 2 from Gary (2001) shows how the plasma f3 value changes with height in the solar atmo- 
sphere. As one can see a region with (3 -C 1 is sandwiched between the photosphere and the upper 
corona, where (3 is about unity or larger. In regions with /3< 1 the magnetic pressure dominates 
over the plasma pressure (and as well over other non-magnetic forces like gravity and the kinematic 
plasma flow pressure) . Here we can neglect in the lowest order all non-magnetic forces and assume 
that the Lorentz force vanishes. This approach is called the force-free field approximation and for 
static configurations it is defined as: 

jxB = 0, (2) 

j = — VxB is the electric current density, (3) 

V B = 0, (4) 

or by inserting Equation (3) into (2): 

(V x B) x B = 0, (5) 

V B = 0. (6) 

Equation (5) can be fulfilled either by: 

V x B = current-free or potential magnetic fields (7) 

or by 

B || V x B force-free fields. (8) 

Current free (potential) fields are the simplest assumption for the coronal magnetic field. The 
line-of-sight (LOS) photosphcric magnetic field which is routinely measured with magnetographs 
are used as boundary conditions to solve the Laplace equation for the scalar potential </>, 

A(j) = (9) 

where the Laplacian operator A is the divergence of the gradient of the scalar field and 

B = -V0. (10) 

When one deals with magnetic fields of a global scale, one usually assumes the so-called "source 
surface" (at about 2.5 solar radii where all field lines become radial): See, e.g., Schatten et al. (1969) 
for details on the potential-field source-surface (PFSS) model. Figure 3 shows such a potential-field 
source-surface model for May 2001 from Wiegelmann and Solanki (2004). 

Potential fields are popular due to their mathematical simplicity and provide a first coarse 
view of the magnetic structure in the solar corona. They cannot, however, be used to model the 
magnetic field in active regions precisely, because they do not contain free magnetic energy to 
drive eruptions. Further, the transverse photospheric magnetic field computed from the potential- 
field assumption usually does not agree with measurements and the resulting potential field lines 
do deviate from coronal loop observations. For example, a comparison of global potential fields 
with TRACE images by Schrijver et al. (2005) and with stereoscopically-reconstructed loops by 
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Figure 3: Global potential field reconstruction. The original figure was published in Wiegelmann 
and Solanki (2004). 



Sandman et al. (2009) showed large deviations between potential magnetic field lines and coronal 
loops. 

The B || V x B condition can be rewritten as 

V x B = qB, (11) 
B-Va = 0, (12) 

where a is called the force-free parameter or force-free function. From the horizontal photospheric 
magnetic field components (B x0 , B y0 ) we can compute the vertical electric current density 

dB y0 dB x0 

and the corresponding distribution of the force-free function a(x, y) in the photosphere 

a(x,y) = no-^-. (14) 

Condition (12) has been derived by taking the divergence of Equation (11) and using the solenoidal 
condition (4). Mathematically Equations (11) and (12) are equivalent to Equations (2) -(4). Pa- 
rameter a can be a function of position, but Equation (12) requires that a be constant along a 
field line. If a is constant everywhere in the volume under consideration, the field is called lin- 
ear force-free field (LFFF), otherwise it is nonlinear force-free field (NLFFF). Equations (11) and 
(12) constitute partial differential equations of mixed elliptic and hyperbolic type. They can be 
solved as a well-posed boundary value problem by prescribing the vertical magnetic field and for 
one polarity the distribution of a at the boundaries. As shown by Bineau (1972) these boundary 
conditions ensure the existence and unique NLFFF solutions at least for small values of a and 
weak nonlinearities. Boulmezaoud and Amari (2000) proved the existence of solutions for a simply 
and multiply connected domain. As pointed out by Aly and Amari (2007) these boundary condi- 
tions disregard part of the observed photospheric vector field: In one polarity only the curl of the 



horizontal field (Equation (13)) is used as the boundary condition, and the horizontal field of the 
other polarity is not used at all. For a general introduction to complex boundary value problems 
with elliptic and hyperbolic equations we refer to Kaiser (2000). 

Please note that high plasma f3 configurations are not necessarily a contradiction to the force- 
free condition (see Neukirch, 2005, for details). If the plasma pressure is constant or the pressure 
gradient is compensated by the gravity force (Vp = — pV^, where p is the mass density and 
the gravity potential of the Sun.) a high-/? configuration can still be consistent with a vanishing 
Lorentz force of the magnetic field. In this sense a low plasma (3 value is a sufficient, but not 
a necessary, criterion for the force-free assumption. In the generic case, however, high plasma (3 
configurations will not be force-free and the approach of the force-free field is limited to the upper 
chromosphere and the corona (up to about 2.5 R s ). 
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2 Linear Force-Free Fields 



Linear force-free fields are characterized by 



V x B = aB, 
VB = 0, 



(15) 
(16) 



where the force-free parameter a is constant. Taking the curl of Equation (15) and using the 
solenoidal condition (16) we derive a vector Helmholtz equation: 



AB + a 2 B = 



(17) 



which can be solved by a separation ansatz, a Green's function method (Chiu and Hilton, 1977) or 
a Fourier method (Alissandrakis, 1981). These methods can also be used to compute a potential 
field by choosing a = 0. 

For computing the solar magnetic field in the corona with the linear force-free model one needs 
only measurements of the LOS photospheric magnetic field. The force-free parameter a is a priori 
unknown and we will discuss later how a can be approximated from observations. Seehafer (1978) 
derived solutions of the linear force-free equations (local Cartesian geometry with (x, y) in the 
photosphere and z is the height from the Sun's surface) in the form: 
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with \ mn = ir 2 (m 2 /Ll + n 2 / ' L 2 y ) and r mn = \/\ mn - a 2 . 

As the boundary condition, the method uses the distribution of B z (x,y) on the photosphere 
z = 0. The coefficients C mn can be obtained by comparing Equation (20) for z — with the 
magnetogram data. In practice Seehafer's (1978) method is used for calculating the linear force- 
free field (or potential field for a = 0) for a given magnetogram (e.g., MDI on SOHO) and a given 
value of a as follows. The observed magnetogram which covers a rectangular region extending 
from to L x in x and to L y in y is artificially extended onto a rectangular region covering — L x 
to L x and — L y to L y by taking an antisymmetric mirror image of the original magnetogram in 
the extended region, i.e., 



B z {-x,y) = -B z (x,y), 
B z {x,-y) = -B z (x,y), 
B z (-x,-y) = B z (x,y) 



(0 < x < L x ,0 < y < L y ). 



This makes the total magnetic flux in the whole extended region to be zero. (Alternatively one 
may pad the extended region with zeros, although in this case the total magnetic flux may be 



7 



non-zero.) The coefficients C mn are derived from this enlarged magnetogram with the help of a 
Fast Fourier Transform. In order for r mn to be real and positive so that solutions (18)-(20) do not 
diverge at infinity, a 2 should not exceed the maximum value for given L x and L y , 



For L x — L y we have L — L x = L y . With this normalization the values of a fall into the range 
-V2tt < a < y/2n. 

2.1 How to obtain the force-free parameter a 

Linear force-free fields require the LOS magnetic field in the photosphere as input and contain a 
free parameter a. One possibility to approximate a is to compute an averaged value of a from the 
measured horizontal photospheric magnetic fields as done, e.g., in Pevtsov et al. (1994), Wheatland 
(1999), Leka and Skumanich (1999), and Hagino and Sakurai (2004), where Hagino and Sakurai 
(2004) calculated an averaged value a = ^2 /ioJ z sign(B z ) / ^2\B Z \. The vertical electric current 

in the photosphere is computed from the horizontal photospheric field as J z = -j- ( — ^Mf- ) ■ 



Such approaches derive best fits of a linear force-free parameter a with the measured horizontal 
photospheric magnetic field. 

Alternative methods use coronal observations to find the optimal value of a. This approach 
usually means that one computes several magnetic field configurations with varying values of a 
in the allowed range and to compute the corresponding magnetic field lines. The field lines are 
then projected onto coronal plasma images. A method developed by Carcedo et al. (2003) is 
shown in Figure 4. In this approach the shape of a number of field lines with different values of 
a, which connect the foot point areas (marked as start and target in Figure 4(e)) are compared 
with a coronal image. For a convenient quantitative comparison the original image shown in 
Figure 4(a) is converted to a coordinate system using the distances along and perpendicular to the 
field line, as shown in Figure 4(b). For a certain number of N points along this uncurled loop the 
perpendicular intensity profile of the emitting plasma is fitted by a Gaussian profile in Figure 4(c) 
and the deviation between field line and loops are measured in Figure 4(d). Finally, the optimal 
linear force-free value of a is obtained by minimizing this deviation with respect to a, as seen in 
Figure 4(f). 

The method of Carcedo et al. (2003) has been developed mainly with the aim of computing 
the optimal a for an individual coronal loop and involves several human steps, e.g., identifying 
an individual loop and its footpoint areas and it is required that the full loop, including both 
footpoints, is visible. This makes it somewhat difficult to apply the method to images with a 
large number of loops and when only parts of the loops are visible. For EUV loops it is also 
often not possible to identify both footpoints. These shortcomings can be overcome by using 
feature recognition techniques, e.g., as developed in Aschwanden et al. (2008a) and Inhester et al. 
(2008) to extract one-dimensional curve-like structures (loops) automatically out of coronal plasma 
images. These identified loops can then be directly compared with the projections of the magnetic 
field lines, e.g., by computing the area spanned between the loop and the field line as defined in 
Wiegelmann et al. (2006b). This method has become popular in particular after the launch of the 
two STEREO spacecraft in October 2006 (Kaiser et at, 2008). The projections of the 3D linear 
force-free magnetic field lines can be compared with images from two vantage viewpoints as done 




Usually a is normalized by the harmonic mean L of L x and L y defined by 
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Figure 4: How to obtain the optimal linear force-free parameter a from coronal observations. [The 
original figures were published as Figures 3, 4, and 5 in Carcedo et al. (2003)]. 
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for example in Feng et al. (2007b, a). This automatic method applied to a number of loops in 
one active region revealed, however, a severe shortcoming of linear force-free field models. The 
optimal linear force-free parameter a varied for different field lines, which is a contradiction to 
the assumption of a linear model. A similar result was obtained by Wiegelmann and Neukirch 
(2002) who tried to fit the loops stereoscopically reconstructed by Aschwanden et al. (1999). On 
the other hand, Marsch et al. (2004) found in their example that one value of a was sufficient to 
fit several coronal loops. Therefore, the fitting procedure tells us also whether an active region can 
be described consistently by a linear force-free field model: Only if the scatter in the optimal a 
values among field lines is small, one has a consistent linear force-free field model which fits coronal 
structures. In the generic case that a changes significantly between field lines, one cannot obtain 
a self-consistent force- free field by a superposition of linear force- free fields, because the resulting 
configurations arc not force-free. As pointed out by Malanushenko et al. (2009) it is possible, 
however, to estimate quantities like twist and loop heights with an error about of 15% and 5%, 
respectively. The price one has to pay is using a model that is not self-consistent. 
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Figure 5: Low and Lou's (1990) analytic nonlinear force-free equilibrium. The original 2D equilib- 
rium is invariant in ip, as shown in panel a. Rotating the 2D-equilibrium and a transformation to 
Cartesian coordinates make this symmetry less obvious (panels b-d), where the equilibrium has 



been rotated by an angle of (p 



8 ' 4 ' 



and respectively. The colour-coding corresponds to the 



vertical magnetic field strength in G (gauss) in the photosphere (z = in the model) and a number 
of arbitrary selected magnetic field lines are shown in yellow. 
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3 Analytic or Semi- Analytic Approaches to Nonlinear Force- 
Free Fields 



Solving the nonlinear force-free equations in full 3-D is extremely difficult. Configurations with 
one or two invariant coordinate(s) are more suitable for an analytic or semi-analytic treatment. 
Solutions in the form of an infinitely long cylinder with axial symmetry are the simplest cases, and 
two best known examples are Lundquist's (1950) solution in terms of Bessel functions ( a = const ant), 
and a solution used by Gold and Hoyle (1960) in their flare model (a ^ constant, all field lines 
have the same pitch in the direction of the axis). Low (1973) considered a ID Cartesian (slab) 
geometry and analyzed slow time evolution of the force-free field with resistive diffusion. 

In Cartesian 2D geometry with one ignorable coordinate in the horizontal (depth) direction, 
one ends up with a second-order partial differential equation, called the Grad-Shafranov equation 
in plasma physics. The force-free Grad-Shafranov equation is a special case of the Grad-Shafranov 
equation for magneto-static equilibria (see Grad and Rubin, 1958), which allow to compute plasma 
equilibria with one ignorable coordinate, e.g. a translational, rotational or helical symmetry. For an 
overview on how the Grad-Shafranov equation can be derived for arbitrary curvilinear coordinates 
with axisymmetry we refer to (Marsh, 1996, section 3.2.). In the cartesian case one finds (see, e.g., 
Sturrock, 1994, section 13.4) 

AA = -\ 2 f(A) (21) 

where the magnetic flux function A depends only on two spatial coordinates and any choice of f(A) 
generates a solution of a magneto-static equilibrium with symmetry. For static equilibria with a 
vanishing plasma pressure gradient the method naturally provides us force-free configurations. A 
popular choice for the generating function is an exponential ansatz, see, e.g. Low(1977), Birn 
et al. (1978), Priest and Milne(1980). The existence of solutions (sometimes multiple, sometimes 
none) and bifurcation of a solution sequence have been extensively investigated (e.g., Birn and 
Schindlcr, 1981). We will consider the Grad-Shafranov equation in spherical polar coordinates in 
the following. 

3.1 Low and Lou's (1990) equilibrium 

As an example we refer to Low and Lou (1990), who solved the Grad-Shafranov equation in 
spherical coordinates (r, 9, ip) for axisymmetric (invariant in cp) nonlinear force-free fields. In this 
case the magnetic field is assumed to be written in the form 

1 / 1 dA dA n \ 
B = — — - — e r - — e e + QeJ, 22 
rant) \r o0 or J 

where A is the flux function, and Q represents the ^-component of the magnetic field B, which 
depends only on A. This ansatz automatically satisfies the solenoidal condition (6), and the force- 
free equation (5) reduces to a Grad-Shafranov equation for the flux function A 

d 2 A l-fi 2 d 2 A ^dQ n 

where \x = cos 9. Low and Lou (1990) looked for solutions in the form 

Q{A) = \A 1+1 ' n {a = ^ ~ A x ' n ) (24) 

with a separation ansatz 

A (r,0) = ^- (25) 
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Here n and A are constants and n is not necessarily an integer; n = 1 and A = corresponds to 
a dipole field. Then Equation (23) reduces to an ordinary differential equation for P(fi), which 
can be solved numerically. Either by specifying n or A, the other is determined as an eigenvalue 
problem (Wolfson, 1995). The solution in 3D space is axisymmetric and has a point source at the 
origin. This symmetry is also visible after a transformation to Cartesian geometry as shown in 
Figure 5(a). The symmetry becomes less obvious, however, when the symmetry axis is rotated with 
respect to the Cartesian coordinate axis; see Figures 5(b)-(d). The resulting configurations are very 
popular for testing numerical algorithms for a 3D NLFFF modeling. For such tests the magnetic 
field vector on the bottom boundary of a computational box is extracted from the semi-analytic 
Low-Lou solution and used as the boundary condition for numerical force-free extrapolations. The 
quality of the reconstructed field is evaluated by quantitative comparison with the exact solution; 
see, e.g., Schrijver et al. (2006). Similarly one can shift the origin of the point source with respect 
to the Sun center and the solution is not symmetric to the Sun's surface and can be used to test 
spherical codes. 

3.2 Titov-Demoulin equilibrium 

Another approach for computing axisymmetric NLFFF solutions has been developed in Titov 
and Demoulin (1999). This model active region contains a current-carrying flux-tube, which is 
imbedded into a potential field. A motivation for such an approach is that solar active regions may 
be thought of as composed of such flux tubes. The method allows to study a sequence of force-free 
configurations through which the flux tube emerges. Figure 6 shows how the equilibrium is built 
up. The model contains a symmetry axis, which is located at a distance d below the photosphere. 
A line current Iq runs along this symmetry axis and creates a circular potential magnetic field. This 
potential field becomes disturbed by a toroidal ring current / with the minor radius a and the major 
radius R, where a <C R is assumed. Two opposite magnetic monopoles of strength q are placed on 
the axis separated by distance L. These monopoles are responsible for the poloidal potential field. 
This field has its field lines overlying the force-free current and stabilizes the otherwise unstable 
configuration. Depending on the choice of parameters one can contain stable or unstable nonlinear 
force-free configurations. The unstable branch of this equilibrium has been used to study the onset 
of coronal mass ejections; see Section 5.5. Stable branches of the Titov-Demoulin equilibrium are 
used as a challenging test for numerical NLFFF extrapolation codes (see, e.g., Wicgelmann et al., 
2006a; Valori et al, 2010). 




Figure 6: Construction of the Titov-Demoulin equilibrium. The original figure was published as 
Figure 2 in Titov and Demoulin (1999). 
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4 Azimuth Ambiguity Removal and Consistency of Field 
Measurements 



4.1 How to derive vector magnetograms? 

NLFFF extrapolations require the photosphcric magnetic field vector as input. Before discussing 
how this vector can be extrapolated into the solar atmosphere, we will address known problems 
regarding the photospheric field measurements. Vector magnetographs are being operated daily at 
NAOJ/Mitaka (Sakurai et at, 1995), NAOC/Huairou (Ai and Hu, 1986), NASA/MSFC (Hagyard 
et ai, 1982), NSO/Kitt Peak (Henney et al, 2006), and U. Hawaii/Mees Observatory (Mickey 
et ai, 1996), among others. The Solar Optical Telescope (SOT; Tsuneta et at, 2008) on the 
Hinode mission has been taking vector magnetograms since 2006. Full-disk vector magnetograms 
are observed routinely since 2010 by the Helioseismic and Magnetic Imager (HMI; Scherrer et al. 
(2012)) onboard the Solar Dynamics Observatory (SDO). Measurements with these vector magne- 
tographs provide us eventually with the magnetic field vector on the photosphere, say B z0 for the 
vertical and B x q and B y o for the horizontal fields. Deriving these quantities from measurements 
is an involved physical process based on the Zeeman and Hanlc effects and the related inversion of 
Stokes profiles (e.g., LaBonte et ai, 1999). Within this work we only outline the main steps and 
refer to del Toro Iniesta and Ruiz Cobo(1996), del Toro Iniesta(2003), and Landi Dcgl'Innocenti 
and Landolfi(2004) for details. Actually measured are polarization degrees across magnetically 
sensitive spectral lines, e.g. the line pair Fe i 6302.5 and 6301.5 A as used on Hinode/SOT (see 
Lites et ai, 2007) or Fe I 6173.3 A as used on SDO/HMI (see Schou et ai, 2012). The accu- 
racy of these measurements depends on the spectral resolution, for example the HMI instruments 
measures at six points in the Fe I 6173.3 A absorption line. In a subsequent step the Stokes 
profiles are inverted to derive the magnetic field strength, its inclination and azimuth. One pos- 
sibility to carry out the inversion (see Lagg et ai, 2004) is to fit the measured Stokes profiles 
with synthetic ones derived from the Unno-Rachkovsky solutions (Unno, 1956; Rachkovsky, 1967). 
Usually one assumes a simple radiative transfer model like the Milne-Eddington atmosphere (see 
e.g. Landi Degl'Innocenti, 1992) in order to derive the analytic Unno-Rachkovsky solutions. The 
line-of-sight component of the field is approximately derived by Bg oc V/I, where V is the circular 
polarization and / the intensity (the so-called weak-field approximation). The error from pho- 
ton noise is approximately SBg oc ^f, where S corresponds to noise in the measured and derived 
quantities. As a rule of thumb, SV/I <~ 10~ 3 and SBg <~ a few gauss (G) in currently operating 
magnetographs. The horizontal field components can be approximately derived from the linear 
polarization Q and U as B 2 oc \J Q 2 + U 2 /I. The error in 5B t is estimated as B t SB t oc Q^Q+ us u 

from which the minimum detectable B t (SB t ~ B t ) is proportional to the square root of the photon 
noise ~ y/ SQ 2 + SU 2 /I w y/SV/I, namely around a few tens of G, one order of magnitude higher 
than SBg. (Although 5B t scales as l/B t and gives much smaller 5B t for stronger B t , one usually 
assumes a conservative error estimate that SB t ~ a few tens of G regardless of the magnitude of 
Bt.) 

Additional complications occur when the observed region is far away from the disk center 
and consequently the line-of-sight and vertical magnetic field components are far apart (see Gary 
and Hagyard, 1990, for details). The inverted horizontal magnetic field components B x0 and B y0 
cannot be uniquely derived, but contain a 180° ambiguity in azimuth, which has to be removed 
before the fields can be extrapolated into the corona. In the following we will discuss this problem 
briefly. For a more detailed review and a comparison and performance check of currently available 
ambiguity-removal routines with synthetic data, see Metcalf et al. (2006). 

To remove the ambiguity from this kind of data, some a priori assumptions regarding the 
structure of the magnetic field vector are necessary, e.g., regarding smoothness. Some methods 
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require also an approximation regarding the 3D magnetic field structure (usually from a potential 
field extrapolation) ; for example to minimize the divergence of magnetic field vector or the angle 
with respect to the potential field. We are mainly interested here in automatic methods, although 
manual methods are also popular, e.g., the AZAM code. If we have in mind, however, the huge 
data stream from SDO/HMI, fully automatic methods are desirable. In the following we will give 
a brief overview on the ambiguity removal techniques and tests with synthetic data. 

4.2 Quantitative comparison of ambiguity removal algorithms 

Metcalf et al. (2006) compared several algorithms and implementations quantitatively with the 
help of two synthetic data sets, a flux-rope simulation by Fan and Gibson (2004) and a multipolar 
constant-a structure computed with the Chiu and Hilton (1977) linear force-free code. The results 
of the different ambiguity removal techniques have been compared with a number of metrics (see 
Table II in Metcalf et at, 2006). For the discussion here we concentrate only on the first test case 
(flux rope) and the area metrics, which simply tells for what fraction of pixels the ambiguity has 
been removed correctly. A value of 1 corresponds to a perfect result and 0.5 to random. The 
result is visualized in Figure 7, where the ambiguity has been removed correctly in black areas. 
Wrong pixels are white. In the following we briefly describe the basic features of these methods 
and provide the performance (fraction of pixels with correctly removed ambiguity) . 

4.3 Ambiguity removal algorithm 

4.3.1 Acute angle method 

The magnetic field in the photosphere is usually not force-free and even not current-free, but an 
often made assumption is that from two possible directions (180° apart) of the observed field B obs , 
the solution with the smaller angle to the potential field (or another suitable reference field) B° 
is the more likely candidate for the true field. Consequently we get for the horizontal/transverse 2 
field components B t the condition 

B° bs • B° > 0. (26) 

This condition is easy to implement and fast in application. In Metcalf et al. (2006) several different 
implementations of the acute angle method are described, which mainly differ by the algorithms 
used to compute the reference field. The different implementations of the acute angle methods got 
a fraction of 0.64-0.75 pixels correct (see Figure 7, panels marked with NJP, YLP, KLP, BBP, 
JLP, and LSPM). 

4.3.2 Improved acute angle methods 

A sophistication of the acute angle method uses linear force-free fields (Wang, 1997; Wang et al., 
2001), where the optimal force-free parameter a is chosen to maximize the integral 

f |B obs • B lff | , N 

S = J \ ohsB J d X dy (27) 

where B lff is the linear force-free reference field. A fraction of 0.87 pixels has been identified 
correctly (see Figure 7 second row, right panel marked with HSO). 

Another approach, dubbed uniform shear method by Moon et al. (2003) uses the acute angle 
method (with a potential field as reference) only as a first approximation and subsequently uses 

2 In the following wc assume observations close to the disk center for simplicity, when the vertical and LOS- 
componcnt arc identical. For observations far away from the disk center one has to resolve first the ambiguity and 
apply coordinate transformations from LOS/transverse to vertical/horizontal fields afterwards. 
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Figure 7: Overview of the performance of different algorithms for removing the 180° azimuth 
ambiguity. The codes have been applied to synthetic data (a flux-rope simulation by Fan and 
Gibson (2004)). In black areas the codes found the correct azimuth and in white areas not. The 
original figure was published as Figure 3 in Metcalf et al. (2006). 
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this result to estimate a uniform shear angle between the observed field and the potential field. 
Then the acute angle method is applied again to resolve the ambiguity, taking into account the 
average shear angle between the observed field and the calculated potential field. A fraction of 
0.83 pixels has been identified correctly. Consequently both methods significantly improve the 
potential-field acute angle method (see Figure 7 third row, center panel marked with USM). 

4.3.3 Magnetic pressure gradient 

The magnetic pressure gradient method (Cuperman et al, 1993) assumes a force-free field and that 
the magnetic pressure B 2 /2 decreases with height. Using the solenoidal and force-free conditions, 
we can compute the vertical magnetic pressure gradient as: 



with any initial choice for the ambiguity of the horizontal magnetic field components (B x ,B y ). 
Different solutions of the ambiguity removal method give the same amplitude, but opposite sign 
for the vertical pressure gradient. If the vertical gradient becomes positive, then the transverse 
field vector is reversed. For the test this method got a fraction of 0.74 pixels correct, which is 
comparable with the potential-field acute angle method (see Figure 7 forth row, left panel marked 
with MS). 

4.3.4 Structure minimization method 

The structure minimization method (Gcorgoulis et al, 2004) is a semi-analytic method which aims 
at eliminating dependencies between pixels. We do not describe the method here, because in the 
test only for a fraction of 0.22 pixels the ambiguity has been removed correctly, which is worse 
than a random result (see Figure 7 third row, right panel marked with MPG). 

4.3.5 Non-potential magnetic field calculation method 

The non-potential magnetic field method developed by Georgoulis (2005) is identical with the acute 
angle method close to the disk center. Away from the disk center the method is more sophisticated 
and uses the fact that the magnetic field can be represented as a combination of a potential field 
and a non-potential part B = B p + B c , where the non-potential part B c is horizontal on the 
boundary and only B c contains electric currents. The method aims at computing a fair a priori 
approximation of the electric current density before the ambiguity removal. With the help of a 
Fourier method the component B c and the corresponding approximate field B arc computed. This 
field is then used as the reference field for an acute angle method. The quality of the reference 
field depends on the accuracy of the a priori assumed electric current density j z . In the original 
implementation by Georgoulis (2005) j z was chosen once a priori and not changed afterwards. 
In an improved implementation (published as part of the comparison paper by Metcalf et al. 
(2006) and implemented by Georgoulis) j z becomes updated in an iterative process. The original 
implementation got 0.70 pixels correct and the improved version 0.90 (see Figure 7 forth row, 
center and right panels marked with NPFC and NPFC2, respectively). So the original method is 
on the same level as the potential-field acute angle method, but the current iteration introduced in 
the updated method gives significantly better results. This method has been used for example to 
resolve the ambiguity of full-disk vector magnetograms from the SOLIS instrument (Hcnney et al., 
2006) at NSO/Kitt Peak. 




(28) 
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4.3.6 Pseudo-current method 



The pseudo-current method developed by Gary and Dcmoulin (1995) uses as the initial step the 
potential-field acute angle method and subsequently applies this result to compute an approx- 
imation for the vertical electric current density. The current density is then approximated by 
a number of local maxima of j z with an analytic expression containing free model parameters, 
which are computed by minimizing a functional of the square of the vertical current density. This 
optimized current density is then used to compute a correction to the potential field. This new 
reference field is then used in the acute angle method to resolve the ambiguity. In the test case this 
method got a fraction of 0.78 of pixels correct, which is only slightly better than the potential- field 
acute angle method (see Figure 7 fifth row, left panel marked with PCM). 

4.3.7 U. Hawai'i iterative method 

This method, originally developed in Canfield et al. (1993) and subsequently improved by a group 
of people at the Institute for Astronomy, U. Hawai'i. As the initial step the acute angle method is 
applied, which is then improved by a constant-a force-free field, where a has to be specified by the 
user (in principle it should also be possible to apply an automatic a-fitting method as discussed 
in Section 4.3.2). Therefore, the result would be similar to the improved acute angle methods, 
but additional two more steps have been introduced for a further improvement. In a subsequent 
step the solution is smoothed (minimizing the angle between neighboring pixels) by starting at 
a location where the field is radial and the ambiguity is obvious, e.g., the umbra of a sunspot. 
Finally also the magnetic field divergence or vertical electric current density is minimized. This 
code includes several parameters, which have to be specified by the user. In the test case the 
code recognized a fraction of 0.97 pixels correctly. So the additional steps beyond the improved 
acute angle method provide another significant improvement and almost the entire region has been 
correctly identified (see Figure 7 fifth row, center panel marked with UHIM). 

4.3.8 Minimum energy methods 

The minimum energy method has been developed by Metcalf (1994). As other sophisticated 
methods it uses the potential-field acute angle method as the initial step. Subsequently a pseudo 
energy, which is defined as a combination of the magnetic field divergence and electric current 
density, is minimized. In the original formulation the energy was defined as E — 5^(|V • B| + |j|), 
which was slightly modified to 



in an updated version. For computing j x , j y , and dB z /dz, a linear force-free model is computed, 
in the same way as described in Section 4.3.7. The method minimizes the functional (29) with 
the help of a simulated annealing method, which is a robust algorithm to find a global minimum. 
In a recent update (published in Metcalf et al., 2006) the (global) linear force-free assumption 
has been relaxed and replaced by local linear force-free assumptions in overlapping parts of the 
magnetogram. The method was dubbed nonlinear minimum energy method, although it does not 
use true NLFF fields (would be too slow) for computing the divergence and electric currents. The 
original linear method got a fraction of 0.98 of pixels correctly and the nonlinear minimum energy 
method even 1.00. Almost all pixels have been correct, except a few on the boundary (see Figure 7 
fifth row, right panel and last row left panel, marked with ME1 and ME2, respectively.) Among 
the fully automatic methods this approach had the best performance on accuracy. A problem 
for practical use of the method was that it is very slow, in particular for the nonlinear version. 
Minimum energy methods are routinely used to resolve the ambiguity in active regions as measured, 
e.g., with SOT on Hinode or HMI on SDO. 




(29) 
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4.4 Summary of automatic methods 

The potential-field acute angle method is easy to implement and fast, but its performance of 
0.64-0.75 is relatively poor. The method is, however, very important as an initial step for more 
sophisticated methods. Using more sophisticated reference fields (linear force-free fields, constant 
shear, non-potential fields) in the acute angle method improves the performance to about 0.83- 
0.90. Linear force-free or similar fields are a better approximation to a suitable reference field, 
but the corresponding assumptions arc not fulfilled in a strict sense, which prevents a higher 
performance. The magnetic pressure gradient and pseudo-current methods are more difficult to 
implement as simple acute angle methods, but do not perform significantly better. A higher 
performance is prevented, because the basic assumptions are usually not fulfilled in the entire 
region. For example, the assumption that the magnetic pressure always decreases with height is 
not fulfilled over bald patches (Titov et al., 1993). The multi-step U. Hawai'i iterative method 
and the minimum energy methods showed the highest performance of > 0.97. The pseudo-current 
method is in principle similar to the better performing minimum energy methods, but due to 
several local minima it is not guaranteed that the method will always find the global minimum. 
Let us remark that Metcalf et al. (2006) introduced more comparison metrics, which, however, 
do not influence the relative rating of the discussed ambiguity algorithms. They also carried out 
another test case using the Chiu and Hilton (1977) linear force-free model, for which most of the 
codes showed an absolutely better performance, but again this does hardly influence the relative 
performance of the different methods. One exception was the improved non-potential magnetic 
field algorithm, which performed with similar excellence as the minimum energy and U. Hawai'i 
iterative methods. Consequently these three methods are all suitable candidates for application to 
data. It is, however, not entirely clear to what extent these methods can be applied to full-disk 
vector magnetograms and what kind of computer resources are required. 

4.4.1 Effects of noise and spatial resolution 

The comparison of ambiguity removal methods started in Metcalf et al. (2006) has been continued 
in Leka et al. (2009). The authors investigated the effects of Poisson photon noise and a limited 
spatial resolution. It was found that most codes can deal well with random noise and the ambiguity 
resolution results are mainly affected locally, but bad solutions (which are locally wrong due to 
noise) do not propagate within the magnetogram. A limited spatial resolution leads to a loss of 
information about the fine structure of the magnetic field and erroneous ambiguity solutions. Both 
photon noise and binning to a lower spatial resolution can lead to artificial vertical currents. The 
combined effect of noise and binning affect the computation of a reference magnetic field used 
in acute angle methods as well as quantities in minimization approaches like the electric current 
density and V • B. Sophisticated methods based on minimization schemes performed again best in 
the comparison of methods and are more suitable to deal with the additional challenges of noise 
and limited resolution. As a consequence of these results Leka et al. (2009) suggested that one 
should use the highest possible resolution for the ambiguity resolution task and if binning of the 
data is necessary, this should be done only after removing the ambiguity. Recently Georgoulis 
(2012) challenged their conclusion that the limited spatial resolution was the cause of the failure of 
ambiguity removal techniques using potential or non-potential reference fields. Georgoulis (2012) 
pointed out that the failure was caused by a non-realistic test-data set and not by the limited 
spatial resolution. This debate has been continued in a reply by Leka et al. (2012). We aim to 
follow the ongoing debate and provide an update on this issue in due time. 
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4.4.2 HAO AZAM method 



This is an interactive tool, which needs human intervention for the ambiguity removal. In the 
test case, which has been implemented and applied by Bruce Lites, all pixels have been identified 
correctly. It is of course difficult to tell about the performance of the method, but only about a 
human and software combination. For some individual or a few active regions the method might 
be appropriate, but not for a large amount of data. 

4.4.3 Ambiguity removal methods using additional observations 

The methods described so far use as input the photospheric magnetic field vector measured at a 
single height in the photosphere. If additional observations/measurements are available they can be 
used for the ambiguity removal. Measurements at different heights in order to solve the ambiguity 
problem have been proposed by Li et al. (1993) and revisited by Li et al. (2007). Knowledge of 
the magnetic field vector at two heights allows us to compute the divergence of the magnetic field 
and the method was dubbed divergence-free method. The method is non- iterative and thus fast. 
Li et al. (2007) applied the method to the same flux-rope simulation by Fan and Gibson (2004) 
as discussed in the examples above, and the method recovered about a fraction of 0.98 pixels 
correctly. The main shortcoming of this method is certainly that it can be applied only if vector 
magnetic field measurements at two heights are available, which is unfortunately not the case for 
most current data sets. 

Martin et al. (2008) developed the so-called chirality method for the ambiguity removal, which 
takes additional observations into account, e.g., Ha, EUV, or X-ray images. Such images are 
used to identify the chirality in solar features like filaments, fibrils, filament channels, or coronal 
loops. Martin et al. (2008) applied the method to different solar features, but to our knowledge the 
method has not been tested with synthetic data, where the true solution of the ambiguity is known. 
Therefore, unfortunately one cannot compare the performance of this method with the algorithms 
described above. It is also now obvious that fully automatic feature recognition techniques to 
identify the chirality from observed images need to be developed. 

After the launch of Solar Orbiter additional vector magnetograms will become available from 
above the ecliptic. Taking these observations from two vantage positions combined is expected to be 
helpful for the ambiguity resolution. If separated by a certain angle, the definition of line-of-sight 
field and transverse field will be very different from both viewpoints. Removing the ambiguity 
should be a straightforward process by applying the transformation to vertical and horizontal 
fields on the photosphere from both viewpoints separately. If the wrong azimuth is chosen, then 
both solutions will be very different and the ambiguity can be removed by simply checking the 
consistency between vertical and horizontal fields from both observations. 

4.5 Derived quantities, electric currents, and a 

The well-known large uncertainties in the horizontal magnetic field component, in particular in 
weak field regions (see Section 4.1), cause large errors when computing the electric current density 
with finite differences via Equation (13). Even more critical is the computation of a with Equa- 
tion (14) in weak field regions and in particular along polarity inversion lines (see e.g., Cuperman 
et al, 1991). The nonlinear force-free coronal magnetic field extrapolation is a boundary value 
problem. As we will see later, some of the NLFFF codes make use of Equation (14) to specify the 
boundary conditions while other methods use the photospheric magnetic field vector more directly 
to extrapolate the field into the corona. 
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4.6 Consistency criteria for force-free boundary conditions 



After Stokes inversion (see Section 4.1) and azimuth ambiguity removal, we derive the photosphcric 
magnetic field vector. Unfortunately there might be a problem, when we want to use these data as 
the boundary condition for NLFFF extrapolations. Due to Metcalf et al. (1995) the solar magnetic 
field is not force-free in the photosphere (finite (3 plasma), but becomes force- free only at about 
400 km above the photosphere. This is also visible in Figure 2 from Gary (2001), which shows 
the distribution of the plasma /? value with height. Consequently the assumption of a force-free 
magnetic field is not necessarily justified in the photosphere. Unless we have information on the 
magnetic flux through the lateral and top boundaries, we have to assume that the photospheric 
magnetic flux is balanced 

/ B z (x,y,0) dxdy = 0, (30) 
Js 

which is usually the case when taking an entire active region as the field of view. In the following 
we review some necessary conditions the magnetic field vector has to fulfill in order to be suitable as 
boundary conditions for NLFFF extrapolations. Molodensky (1969, 1974) and Aly (1989) defined 
several integral relations, which are related to two moments of the magnetic stress tensor. 

1. The first moment corresponds to the net magnetic force, which has to vanish on the boundary: 

/ B x B z dxdy= / B y B z dxdy = 0, (31) 
Js Js 

I {B 2 X + B 2 ) dx dy = [ B 2 z dxdy. (32) 
Js Js 

2. The second moment corresponds to a vanishing torque on the boundary: 

x{Bl + B 2 )dxdy = xB 2 dxdy, (33) 
Js Js 

[ y(B 2 x + B 2 y )dxdy = [ yB 2 z dxdy, (34) 
Js Js 

/ yB x B z dxdy = / xB y B z dxdy. (35) 
Js Js 



The total energy of a force-free configuration can be estimated directly from boundary conditions 
with the help of the virial theorem (see, e.g., Aly, 1989, for a derivation of this formula) 

E tot = — I (x B x + y B v ) B z dx dy. (36) 
Mo Js 

For Equation (36) to be applicable, the boundary conditions must be compatible with the force-free 
assumption. If the integral relations (31) -(35) are not fulfilled then the data are not consistent 
with the assumption of a force-free field. A principal way to avoid this problem would be to measure 
the magnetic field vector in the low-/? chromosphere, but unfortunately such measurements are not 
routinely available. We have therefore to rely on photospheric measurements and apply some 
procedure, dubbed 'preprocessing', in order to derive suitable boundary conditions for NLFFF 
extrapolations. As pointed out by Aly (1989) the condition that a is constant on magnetic field 
lines (12) leads to the integral relation 

f f(a)B n -dA= f f(a)B n -dA, (37) 
Js+ Js- 
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where S + and 5_ correspond to areas with positive and negative B z in the photosphere, respec- 
tively, and / is an arbitrary function. Condition (37) is referred to as differential flux-balance 
condition as it generalizes the usual flux-balance condition (30). As the connectivity of magnetic 
field lines (magnetic positive and negative regions on the boundary connected by field lines) is 
a priori unknown, relation (37) is usually only evaluated after a 3D force-free model has been 
computed. 

4.7 Preprocessing 

Wiegclmann et al. (2006b) developed a numerical algorithm in order to use the integral rela- 
tions (31) -(35) to derive suitable NLFFF boundary conditions from photosphcric measurements. 
To do so, we define the functional: 

iprop = fJ>iLi + H2L2 + H3L3 + M4-^4, (38) 
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The first and second terms (L\, L2) arc quadratic forms of the force and torque balance conditions, 
respectively. The L 3 term measures the difference between the measured and prcprocessed data. 
L 4 controls the smoothing, which is useful for the application of the data to finite-difference nu- 
merical code and also because the chromospheric low-/? field is smoother than in the photosphere. 
The aim is to minimize L prcp so that all terms L n are made small simultaneously. The optimal 
parameter sets \i n have to be specified for each instrument separately. The resulting magnetic 
field vector is then used to prescribe the boundary conditions for NLFFF extrapolations. In an 
alternative approach Fuhrmann et al. (2007) applied a simulated annealing method to minimize 
the functional. Furthermore they removed the L3 term in favor of a different smoothing term 
L4, which uses the median value in a small window around each pixel for smoothing. The pre- 
processing routine has been extended in Wiegelmann et al. (2008) by including chromospheric 
measurements, e.g., by minimizing additionally the angle between the horizontal magnetic field 
and chromospheric Ha fibrils. In principle one could add additional terms to include more direct 
chromospheric observations, e.g., line-of-sight measurements of the magnetic field in higher regions 
as provided by SOLIS. In principle it should be possible to combine methods for ambiguity removal 
and preprocessing in one code, in particular for ambiguity codes which also minimize a functional 
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like the Metcalf (1994) minimum energy method. A mathematical difficulty for such a combination 
is, however, that the preprocessing routines use continuous values, but the ambiguity algorithms 
use only two discrete states at each pixel. Preprocessing minimizes the integral relations (31-35) 
and the value of these integrals reduces usually by orders of magnitudes during the preprocessing 
procedure. These integral relation are, however, only necessary and not sufficient conditions for 
force- free consistent boundary conditions, and preprocessing does not make use of condition (37). 
Including this condition is not straight forward as one needs to know the magnetic field line con- 
nectivity, which is only available after the force- free configuration has been computed in 3D. An 
alternative approach for deriving force-free consistent boundary conditions is to allow changes of 
the boundary values (in particular the horizontal field) during the force-free reconstruction itself, 
e.g., as recently employed by Wheatland and Rcgnicr (2009), Amari and Aly (2010), and Wiegel- 
mann and Inhester (2010). The numerical implementation of these approaches does necessarily 
depend on the corresponding force-free extrapolation codes and we refer to Sections 6.2 and 6.4 
for details. 

5 Nonlinear Force-free Fields in 3D 

In the following section we briefly discuss some general properties of force-free fields, which are 
relevant for solar physics, like the magnetic helicity, estimations of the minimum and maximum 
energy a force- free field can have for certain boundary conditions and investigations of the stability. 
Such properties are assumed to play an important role for solar eruptions. The Sun and the solar 
corona are of course three-dimensional and for any application to observed data, configurations 
based on symmetry assumptions (as used in Section 3) are usually not applicable. The numerical 
treatment of nonlinear problems, in particular in 3D, is significantly more difficult than linear 
ones. Linearized equations arc often an over-simplification which does not allow the appropriate 
treatment of physical phenomena. This is also true for force-free coronal magnetic fields and has 
been demonstrated by comparing linear force- free configurations (including potential fields, where 
the linear force- free parameter a is zero). 

Computations of the photospheric a distribution from measured vector magnetograms by Equa- 
tion (14) show that a is a function of space (see, e.g., Pevtsov et at, 1994; Regnier et at, 2002; 
DeRosa et ai, 2009). Complementary to this direct observational evidence that nonlinear effects 
are important, there are also theoretical arguments. Linear models are too simple to estimate 
the free magnetic energy. Potential fields correspond to the minimum energy configuration for a 
given magnetic flux distribution on the boundary. Linear force-free fields contain an unbounded 
magnetic energy in an open half-space above the photosphere (Scehafer, 1978), because the govern- 
ing equation in this case is the Hclmholtz (wave) equation (Equation (17)) whose solution decays 
slowly toward infinity. Consequently both approaches are not suitable for the estimation of the 
magnetic energy, in particular not an estimation of the free energy a configuration has in excess of 
a potential field. 

5.1 Magnetic helicity 

Magnetic helicity is a quantity closely related to a property of the force-free field (Woltjer, 1958), 
and is defined by 



where B = V x A and A is the vector potential. When B is given, A is not unique and a gradient 
of any scalar function can be added without changing B. Such gauge freedom does not affect 
the value of H m if the volume V is bounded by a magnetic surface (i.e., no field lines go through 




(43) 
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the surface). Figure 8 shows simple torus configurations and their magnetic helicities. As can 
be guessed from the figures, magnetic helicity is a topological quantity describing how the field 
lines are twisted or mutually linked, and is conserved when resistive diffusion of magnetic field 
is negligible. In the case of the solar corona, the bottom boundary (the photosphere) is not a 
magnetic surface, and field lines go through it. Even under such conditions, an alternative form for 
the magnetic helicity which does not depend on the gauge of A can be defined (Berger and Field, 
1984; Finn and Antonsen, 1985). On the Sun one finds the hemispheric helicity sign rule (see, 
e.g. Pevtsov et at, 1995; Wang and Zhang, 2010, and references therein). For various features like 
active regions, filaments, coronal loops and interplanetary magnetic clouds the helicity is negative 
in the northern and positive in the southern hemisphere. 






H m =0 



H m =ro 2 



H m =±20,<t„ 



Figure 8: Magnetic helicity of field lines in torus configuration: untwisted (left), twisted by T turns 
(middle), and two untwisted but intersecting tori (right). $ stands for the total magnetic flux. 



5.2 Energy principles 

Energy principles leading to various magnetic fields (potential fields, linear force-free fields, and 
nonlinear force-free fields) were summarized in Sakurai (1989). For a given distribution of magnetic 
flux (B z ) on the boundary, 

(a) a potential field is the state of minimum energy. 

(b) If the magnetic energy is minimized with an additional condition of a fixed value of H m , one 
obtains a linear force-free field. The value of constant a should be an implicit function of 
i/ m . The obtained solution may or may not be a minimum of energy; in the latter case the 
solution is dynamically unstable. 

(c) If the magnetic energy is minimized by specifying the connectivity of all the field lines, one 
obtains a nonlinear force-free field. The solution may or may not be dynamically stable. 

Item (c) is more explicitly shown by introducing the so-called Euler potentials (u, v) for the mag- 
netic field (Stern, 1970), 

B = Vu x Vv. (44) 

This representation satisfies V • B = 0. Since B • Vu = B ■ Vw = 0, it and v are constant along 
the field line. The values of u and v on the boundary can be set so that B z matches the given 
boundary condition. If the magnetic energy is minimized with the values of u and v specified on 
the boundary, one obtains Equation (5) for a general (nonlinear) force-free field. 

By the construction of the energy principles, the energy of (b) or (c) is always larger than that 
of the potential field (a). If the values of u and v are so chosen (there is enough freedom) that 
the value of H m is the same in cases (b) and (c), then the energy of nonlinear force- free fields (c) 
is larger than that of the linear force- free field (b). Therefore, we have seen that magnetic energy 
increases as one goes from a potential field to a linear force-free field, and further to a nonlinear 
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force-free field. Suppose there are field lines with enhanced values of a (carrying electric currents 
stronger than the surroundings) . By some instability (or magnetic reconnection) , the excess energy 
may be released and the twist in this part of the volume may diminish. However, in such rapid 
energy release processes, the magnetic helicity over the whole volume tends to be conserved (Berger, 
1984). Namely local twists represented by spatially- varying a only propagate out from the region 
and are homogenized, but do not disappear. Because of energy principle (b), the end state of 
such relaxation will be a linear force-free field. This theory (Taylor relaxation; Taylor, 1974, 1986) 
explains the commonly-observed behavior of laboratory plasmas to relax toward linear force-free 
fields. On the Sun this behaviour is not observed, however. A possible explanation could be that 
since we observe spatially-varying a on the Sun, relaxation to linear force-free fields only takes 
place at limited occasions (e.g., in a flare) and over a limited volume which magnetic reconnection 
(or other processes) can propagate and homogenize the twist. 

5.3 Maximum energy 

There is in particular a large interest on force-free configurations for a given vertical magnetic 
field B n on the lower boundary and in which range the energy content of these configurations 
can be. For such theoretical investigations, one usually assumes a so-called star-shaped volume, 
like the exterior of a spherical shell and the coronal magnetic field is unbounded but has a finite 
magnetic energy. (Numerical computations, on the other hand, are mainly carried out in finite 
computational volumes, like a 3D-box in Cartesian geometry.) It is not the aim of this review 
to follow the involved mathematical derivation, which the interested reader finds in Aly (1984). 
As we saw above, the minimum energy state is reached for a potential field. On the other hand, 
one is also interested in the maximum energy a force-free configuration can obtain for the same 
boundary conditions B n . This problem has been addressed in the so-called Aly-Sturrock conjecture 
(Aly, 1984, 1991; Sturrock, 1991). The conjecture says that the maximum magnetic energy is 
obtained if all magnetic field lines are open (have one footpoint in the lower boundary and reach 
to infinity). This result implies that any non-open force-free equilibrium (which contains electric 
currents parallel to closed magnetic field lines, e.g., created by stressing closed potential field lines) 
contains an energy which is higher than the potential field, but lower than the open field. As 
pointed out by Aly (1991) these results imply that the maximum energy which can be released 
from an active region, say in a flare or coronal mass ejection (CME), is the difference between 
the energy of an open field and a potential field. While a flare requires free magnetic energy, 
the Aly-Sturrock conjecture does also have the consequence that it would be impossible that 
all field lines become open directly after a flare, because opening the field lines costs energy. 
This is in some way a contradiction to observations of CMEs, where a closed magnetic structure 
opens during the eruption. Choe and Cheng (2002) constructed force-free equilibria containing 
tangential discontinuities in multiple flux systems, which can be generated by footpoint motions 
from an initial potential field. These configurations contain energy exceeding the open field, a 
violation of the Aly-Sturrock conjecture, and would release energy by opening all field lines. Due 
to the tangential discontinuities, these configurations contain thin current sheets, which can develop 
micro-instabilities to convert magnetic energy into other energy forms (kinetic and thermal energy) 
by resistive processes like magnetic reconnection. It is not clear (Aly and Amari, 2007), however, 
which conditions are exactly necessary to derive force-free fields with energies above the open field: 
Is it necessary that the multiple flux-tubes are separated by non-magnetic regions like in Choe and 
Cheng (2002)? Or would it be sufficient that the field in this region is much weaker than in the 
flux tubes but remains finite? (See Sakurai, 2007, for a related discussion.) 
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5.4 Stability of force-free fields 



In principle the MHD stability criteria can also be applied to force-free equilibria. Typical ap- 
proaches (see the book by Priest, 1982) to investigate the stability of ideal MHD equilibria (which 
correspond to the assumption of infinite electrical conductivity) are normal mode analysis and an 
energy criterion. The basic question is how a small disturbance to the equilibrium evolves. Analytic 
methods typically linearize the problem around an equilibrium state, which leads to the so-called 
linear stability analysis. One has to keep in mind, however, that a linearly-stable configuration 
might well be nonlinearly unstable. The nonlinear stability of a system is usually investigated 
numerically with the help of time dependent simulations, e.g., with an MHD code (see also Sec- 
tion 5.5 for an application to NLFFF equilibria). In the following we concentrate on linear stability 
investigations by using an energy criterion. 

For a force-free configuration the energy is given by 

Wo = [ ^dV (45) 
J 2^o 

where the subscript corresponds to the equilibrium state. This equilibrium becomes disturbed 
by an displacement £(ro,f) in the form B = Bo + Bi with Bi = Vo x (£ x Bo). This form 
of the magnetic field displacement has its origin from the linearized induction equation = 
V x (vi x Bo), where the velocity field has been replaced by the displacement £. The MHD energy 
principle (Bernstein et at, 1958) reduces for force- free fields to (Molodcnsky, 1974): 

W = — [ [(V x (£ x B)) 2 - (V x (£ x B)) ■ (£ x (V x B))l dV. (46) 
2^0 Jv L " J 

A configuration is stable if W > 0, unstable for W < and marginally stable for W — 0. For 
force-free fields and using the perturbed vector potential Ai = £ x B, Equation (46) can be written 
as: 

2/xo J v 



(V x AO 2 -aAi-VxAi dV. (47) 



From Equation (47) it is obvious that the potential field with a — is stable. If we approximate 
|V x Ai| <~ \Ai\/£ with a typical length scale I of the system, the first term may remain larger 
than the second term (i.e., stability) in Equation (47) if 

\a\<\ji. (48) 

This means that the scale of twist in the system, l/a, should be larger than the system size I 
for it to be stable. This criterion is known as Shafranov's limit in plasma physics. More precise 
criteria for stability can be obtained for specific geometries. For example the case of cylindrical 
linear force-free field (Lundquist's field) was studied by Goedbloed and Hagebeuk (1972). 



5.5 Numerical stability investigations 

Tdrok and Kliem (2005) investigated the stability of the nonlinear force-free Titov-Demoulin equi- 
librium numerically with the help of a time-dependent MHD code. Figure 9 shows snapshots 
from MHD simulation starting from in an unstable branch of the Titov-Demoulin equilibrium in 
comparison with a solar eruption observed with TRACE. The simulation shows a very good agree- 
ment with the observed eruptions and indicates that a helical kink instability can trigger coronal 
eruptions. Dependent on particular parameters of the original Titov-Demoulin equilibrium the 
eruption remains confined or leads to a coronal mass ejection (see Torok and Kliem, 2005, for 
details) . 
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Figure 9: Numerical simulations starting from an unstable branch of the Titov-Demoulin equilib- 
rium in comparison with TRACE observations of an eruption. The original figure was published 
as Figure 1 in Torok and Kliem (2005). 
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6 Numerical Methods for Nonlinear Force-free Fields 



In the following we review five different approaches for the computation of nonlinear force-free 
coronal magnetic fields. The aim of all codes is to extrapolate photospheric vector field measure- 
ments into the corona, but the way how the measurements are used is different. MHD relaxation 
and optimization methods prescribe the three components of the magnetic field vector on the 
bottom boundary. Grad-Rubin methods use the vertical magnetic field and the vertical electric 
current density (or a-distribution) as boundary condition. The upward integration method and 
the boundary-element method require a combination of both boundary conditions. In the follow- 
ing we will briefly discuss the main features of these five methods. Grad-Rubin, MHD relaxation, 
and optimization methods require first the computation of a potential field; then the appropriate 
boundary conditions are specified and eventually one iterates numerically for a solution of the 
NLFFF equations. Upward integration and boundary-element methods do not require first the 
computation of a potential field, but solve the NLFFF equations more directly. Both methods 
have, however, some shortcomings as explained later. Often one is interested anyway to get also 
the potential field, e.g., to derive the energy the NLFFF field has in excess of the potential field. 
A more detailed review on the mathematical and computational implementations and recent code 
updates can be found in Wiegelmann (2008). 

6.1 Upward integration method 

This straightforward method was proposed by Nakagawa (1974) and it has been first computa- 
tionally implemented by Wu et al. (1985, 1990). The basic idea of this method is to reformulate 
Equations (2) - (4) and extrapolate the magnetic field vector into the solar corona. The method is 
not iterative and extrapolates the magnetic field directly upward, starting from the bottom layer, 
where the field is measured. From Bo (a;, y, 0) one computes the z-componcnt of the electric current 
MoJzO by Equation (13) and the corresponding a-distribution with Equation (14). Then the x and 
y-components of the electric current are calculated by Equation (11): 

Moixo = a Q B x0 , (49) 
Mojyo = a a B y0 . (50) 

Finally, we get the z-derivatives of the magnetic field vector with Equations (3) and (4) as 

dB x0 . dB z0 , . 

= VOjyO + , (51) 



dz ^ Jy " dx 

dB y0 dB z0 



Mojxo , (52) 



dz dy 

(53) 



9B z q dB X Q dByo 



dz dx dy 

A numerical integration provides the magnetic field vector at the level z + dz. These steps are re- 
peated in order to integrate the equations upwards in z. Naively one would assume to derive finally 
the 3D magnetic fields in the corona, which is indeed the idea of this method. The main problem 
is that this simple straightforward approach does not work because the method is mathematically 
ill-posed and the algorithm is unstable (see, e.g., Cuperman et al., 1990 and Amari et al., 1997 for 
details). As a result of this numerical instability one finds an exponential growth of the magnetic 
field w ith increasing height. The reason for this is that the method transports information only 
from the photosphere upwards. Other boundary conditions, e.g., at an upper boundary, either at 
a finite height or at infinity cannot be taken into account. Several attempts have been made to 
stabilize the algorithm, e.g., by smoothing and reformulating the problem with smooth analytic 
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functions (e.g., Cuperman et al., 1991; Demoulin and Priest, 1992; Song et al, 2006). Smoothing 
does help somewhat to diminish the effect of growing modes, because the shortest spatial modes 
are the fastest growing ones. To our knowledge the upward integration method has not been com- 
pared in detail with other NLFFF codes and it is therefore hard to evaluate the performance of 
this method. 



6.2 Grad Rubin method 

The Grad-Rubin method has been originally proposed (but not numerically implemented) by Grad 
and Rubin (1958) for the application to fusion plasma. The first numerical application to coronal 
magnetic fields was carried out by Sakurai (1981). The original Grad-Rubin approach uses the a- 
distribution on one polarity and the initial potential magnetic field to calculate the electric current 
density with Equation (12) and to update the new magnetic field B from the Biot-Savart equation 
(11). This scheme is repeated iteratively until a stationary state is reached, where the magnetic 
field does not change anymore. Amari et al. (1997, 1999) implemented the Grad-Rubin method 
on a finite difference grid and decomposes Equations (2) - (4) into a hyperbolic part for evolving a 
along the magnetic field lines and an elliptic one to update the magnetic field from Ampere's law: 

BW.VaW = 0, (54) 
a {k) \ s ± = ao±- (55) 

This evolves a from one polarity on the boundary along the magnetic field lines into the volume 
above. The value of «o± is given either in the positive or negative polarity: 

V x B( fe+1 ) = a< fc >B< fc > (56) 

V-B( fe+1 ) = 0, (57) 

£( fc+1 V = B z0 , (58) 

lim |B< fe+1 )| = 0. (59) 

|r|— >oo 

An advantage from a mathematical point of view is that the Grad-Rubin approach solves the 
nonlinear force-free equations as a well-posed boundary value problem. As shown by Bineau 
(1972) the Grad-Rubin-type boundary conditions, the vertical magnetic field and for one polarity 
the distribution of a, ensure the existence and unique NLFFF solutions at least for small values 
of a and weak nonlinearities. See Amari et al. (1997, 2006) for more details on the mathematical 
aspect of this approach. The largest-allowed current and the corresponding maximum values 
of a for which one can expect convergence of the Grad-Rubin approach have been studied in 
Inhester and Wiegelmann (2006). Starting from an initial potential field the NLFFF equations are 
solved iteratively in the form of Equations (11) -(12). The horizontal component of the measured 
magnetic field is then used to compute the distribution of a on the boundary using Equation (14). 
While a is computed this way on the entire lower boundary, the Grad-Rubin method requires 
only the prescription of a for one polarity. For measured data which contain noise, measurement 
errors, finite forces, and other inconsistencies the two solutions can be different: However, see for 
example the extrapolations from Hinode data carried out in DeRosa et al. (2009). While both 
solutions are based on well-posed mathematical problems, they are not necessary consistent with 
the observations on the entire lower boundary. One can check the consistency of the a-distribution 
on both polarities with Equation (37). 

As a further step to derive one unique solution the Grad-Rubin approach has been extended 
by Wheatland and Regnier (2009) and Amari and Aly (2010) by using these two different solutions 
(from different polarities) to correct the a-distribution on the boundaries and to find finally one 
consistent solution by an outer iterative loop, which changes the a-distribution on the boundary. 
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An advantage in this approach is that one can specify where the a-distribution, as computed by 
Equation (14), is trustworthy (usually in strong field regions with a low measurement error in the 
transverse field) and where not (in weak field regions). This outer iterative loop, which aims at 
finding a consistent distribution of a on both polarities, allows also to specify where the initial 
distribution of a is trustworthy. 

6.3 MHD relaxation method 

MHD relaxation method means that a reduced set of time-dependent MHD equations is used to 
compute stationary equilibria: 

vw = (V x B) x B, (60) 
e + vxB = 0, (61) 

— = -Vxe, (62) 
V B = 0. (63) 

Here v is a fictitious viscosity, v the fluid velocity and e the electric field. For general MHD 
equilibria the approach was proposed by Chodura and Schlutcr (1981). Applications to force-free 
coronal magnetic fields can be found in Mikic and McClymont (1994), Roumcliotis (1996), and 
McClymont et al. (1997). In principle any time-dependent MHD code can be used for this aim. 
The first NLFFF implementation of this methods used the code developed by Mikic et al. (1988). 
MHD relaxation means that an initial non-equilibrium state is relaxed towards a stationary state, 
here NLFFF. The initial non-equilibrium state is often a potential field in the 3D-box, where the 
bottom boundary field has been replaced by the measurements. This leads to large deviations from 
the equilibrium close to this boundary. As a consequence one finds a finite plasma flow velocity v in 
Equation (60) because all non-magnetic forces accumulate in the velocity field. This velocity field 
is reduced during the relaxation process and the force-free field equations are obviously fulfilled 
when the left-hand side of Equation (60) vanishes. The viscosity v is usually chosen as 

z/=-|B| 2 (64) 
A 4 

with \i = constant. By combining Equations (60), (61), (62), and (64) one gets a relaxation process 
for the magnetic field 

(9B 

— =/jF mhd , (65) 

F MHD = Vx( " V * B ^ B ' XB ). (66, 
For details regarding a currently- used implementation of this approach see Valori et al. (2005). 

6.4 Optimization approach 

The optimization approach as proposed in Wheatland et al. (2000) is closely related to the MHD 
relaxation approach. It shares with this method that a similar initial non-equilibrium state is 
iterated towards a NLFFF equilibrium. It solves a similar iterative equation as Equation (65) 

£-,F, (67, 
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but F has additional terms, as explained below. The force-free and solenoidal conditions are solved 
by minimizing the functional 

L= [ [B- 2 |(VxB)xB| 2 + |V-B| 2 ] dV. (68) 
Jv 

If the minimum of this functional at L = is attained then the NLFFF equations (2) - (4) are 
fulfilled. The functional is minimized by taking the functional derivative of Equation (68) with 
respect to an iteration parameter t: 



1 dL f dB ^ f dB „ jn 

- — = - -— ■ F dV — / — • G dS, (69) 

2 dt J v dt Js dt y ' 



V x 



[(V x B) x B] x B 



+ -V x 



B 2 

((V -B) B) x B 

-n x (v x b) - v(n • b) 

+ ft(V-B) + ft 2 B}, (70) 

n = B- 2 [(V x B) x B - (V • B) B] . (71) 

For vanishing surface terms the functional L decreases monotonically if the magnetic field is iterated 
by 

f=,F. (72) 

The first term in Equation (70) is identical with Fmhs as defined in Equation (66). 

A principal problem with the optimization and the MHD-relaxation approaches is that using 
the full magnetic field vector on the lower boundary does not guarantee the existence of a force-free 
configuration (see the consistency criteria in Section 4.6). Consequently, if fed with inconsistent 
boundary data, the codes cannot find a force-free configuration, but a finite residual Lorentz force 
and/or a finite divergence of the field remains in the 3D equilibrium. A way around this problem is 
to preprocess the measured photospheric data, as explained in Section 4.7. An alternative approach 
is that one allows deviations of the measured horizontal field vector and the corresponding field 
vector on the lower boundary of the computational box during the minimization of the functional 
(68). Wiegelmann and Inhester (2010) extended this functional by another term 

v [ (B - B obs ) W (B B obs ) dS (73) 
Js 

where v is a free parameter and the matrix W contains information how reliable the data (mainly 
measurements of the horizontal photospheric field) are. With this approach inconsistencies in the 
measurement lead to a solution compatible with physical requirements (vanishing Lorentz force and 
divergence), leaving differences between B Q b s and the bottom boundary field B in regions where W 
is low (and the measurement error high). Consequently this approach takes measurement errors, 
missing data, and data inconsistencies into account. Further tests are necessary to investigate 
whether this approach or preprocessing, or a combination of both, is the most effective way to deal 
with noisy and inconsistent photospheric field measurements. This approach, as well as a variant 
of the Grad-Rubin method, have been developed in response to a joint study by DeRosa et al. 
(2009), where one of the main findings w as that force-free extrapolation codes should be able to 
incorporate measurement inconsistencies (see also Section 6.6). 
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6.5 Boundary-element methods 



The boundary-element method was developed by Yan and Sakurai (2000) and requires the magnetic 
field vector and the a-distribution on the boundary as input. The NLFFF equations relate the 
magnetic field values on the boundary with those in the volume: 

■a* (74) 

with Cj = 1 for points in the volume and q = 1/2 for boundary points and B is the magnetic field 
vector on the boundary, where 

Y = diag ( C0S ( Aa;r ) C0S (V) C0S ( A * r A / 75 ) 
and Aj (i = x, y, z) are implicitly computed with integrals over the 3D volume, 

/ y.^-a^-lVaxB,)^ = 0. (76) 
Jv 

The boundary-element method is slow for computing the NLFFF in a 3D domain. Rudenko and 
Myshyakov (2009) raised questions on this method. 



6.6 Comparison of methods and the NLFFF consortium 

Since 2004 a group of scientists chaired by Karcl Schrijver compare, evaluate, and improve methods 
for the nonlinear force-free computation of coronal magnetic fields and related topics. The test 
cases are available at http://www.lmsal.com/-~derosa/for_nlfrf/ So far six workshops have been 
organized and the consortium published four joint publications: 

1. Schrijver et al. (2006) performed blind tests on analytical force-free field models with various 
boundary conditions to show that in general the NLFFF algorithms perform best where the 
magnetic field and the electrical currents are strongest, but they are also very sensitive to 
the specified boundary conditions. Nevertheless, it was shown that the optimization method 
as proposed by Wheatland et al. (2000) and as implemented by Wiegelmann (2004) was the 
fastest-converging and best-performing one for this analytical test case. 

2. Metcalf et al. (2008) tested the performance of the NLFFF algorithms applied to a solar-like 
reference model including realistic photospheric Lorentz forces and a complex magnetic field 
structure. All the codes were able to recover the presence of a weakly twisted, helical flux 
rope. Due to the sensitivity to the numerical details, however, they were less accurate in 
reproducing the field connectivity and magnetic energy when applied to the preprocessed, 
more force-free, chromospheric-likc boundary conditions. When applied to the forced, not 
preprocessed photospheric data the codes did not perform successfully, indicating that the 
consistency of the used boundary conditions is crucial for the success of the magnetic field 
extrapolations. It also showed that the magnetic field connection between the photosphere, 
chromosphere, and lower corona needs to be additionally precisely modeled. 

3. Schrijver et al. (2008) used four different codes and a variety of boundary conditions to 
compute 14 NLFFF models based on Hinode/SOT-SP 3 data of an active region around 
the time of a powerful flare. When applied to this real solar data, the models produced 
a wide variety of magnetic field geometries, energy contents, and force-freeness. Force- free 

3 Solar Optical Telescope Spectro-Polarimeter 
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consistency criteria, like the alignment of electric currents with magnetic field lines, have 
been best fulfilled for computations with the Grad-Rubin approach. It was concluded that 
strong electrical currents in the form of an ensemble of thin strands emerge together with 
magnetic flux preceding the flare. The global patterns of magnetic fields are compatible with 
a large-scale twisted flux rope topology, and they carry energy which is large enough to power 
the flare and its associated CME. 

4. DcRosa et al. (2009) found that various NLFFF models differ remarkably in the field line 
configuration and produce different estimates of the free magnetic energy when applied to 
Hinodc/SOT-SP data. This problem was recognized already in the first application to Hinode 
data in Schrijvcr et al. (2008) and it has been worked out that a small field-of-view vector 
magnetogram, which docs not contain an entire active region and its surroundings, does 
not provide the necessary magnetic connectivity for successful NLFFF extrapolations. As 
visible in Figure 10 the stereoscopically-reconstructed loops by Aschwanden et al. (2008b) do 
not agree well with the NLFFF models. Unfortunately, the FOV of Hinode covered only a 
small fraction (about 10%) of area spanned by loops reconstructed from STEREO/SECCHI 
images. The quantitative comparison was unsatisfactory and NLFFF models have not been 
better as potential fields here. In other studies NLFFF methods have shown to be superior 
to potential and linear force- free extrapolations (Wiegelmann et at, 2005). NLFF field lines 
showed in particular excellent agreement with the observed loops, when both footpoints are 
within the FOV of the vector magnetogram and sufficiently far away from the boundaries. 

When presented with complete and consistent boundary conditions, NLFFF algorithms gener- 
ally succeed in reproducing the test fields. However, for a well-observed dataset (a Hinode/SOT-SP 
vector-magnetogram embedded in MDI data) the NLFFF algorithms did not yield consistent solu- 
tions. From this study we conclude that one should not rely on a model-field geometry or energy 
estimates unless they match coronal observations. It was concluded that successful application to 
real solar data likely requires at least: 

1. Large model volumes with high resolution that accommodate most of the field- line connec- 
tivity within a region and to its surroundings. 

2. Accommodation of measurement uncertainties (in particular in the transverse field compo- 
nent) in the lower boundary condition. 

3. 'Preprocessing' of the lower-boundary vector field that approximates the physics of the 
photosphere-to-chromosphere interface as it transforms the observed, forced, photospheric 
field to a realistic approximation of the chromospheric, nearly-force-free, field. 

4. The extrapolated coronal magnetic field lines should be compared and verified by coronal 
observations. 

In the meantime some work has been done in reply to these conclusions. New implementations 
of the Grad-Rubin and optimization methods do accommodate the measurement errors; see Sec- 
tions 6.2 and 6.4 for an overview and Wheatland and Regnier (2009), Wiegelmann and Inhester 
(2010), and Amari and Aly (2010) for the corresponding original publications. On the instrumen- 
tation side SDO/HMI provides us with full-disk measurements of the photospheric magnetic field 
vector, which should allow us to find suitable large model volumes. The first vector magnetograms 
from SDO/HMI have been released at the end of 2011 and currently research on using them for 
force-free extrapolations is ongoing. 
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Figure 10: A series of coaligned images of AR 10953. Blue lines are stereoscopically-reconstructed 
loops (Aschwanden et al., 2008b). Red lines are extrapolated nonlinear force-free field lines from 
Hinode/SOT with MDI data outside the Hinode FOV (the dotted line). The original figure was 
published as Figure 1 in DeRosa et al. (2009). 

6.7 Application of nonlinear force-free codes 

Despite the difficulties outlined in Section 6.6 NLFFF-codes have been used to study active regions 
in various situations. Several studies deal with the energy content of the coronal magnetic field in 
active regions. Bleybel et al. (2002) studied the energy budget of AR 7912 before and after a flare 
on 1995 October 14 with a Grad-Rubin method and found that the magnetic energy decreased 
during the flare. The magnetic field lines computed from the nonlinear force-free model seem to 
be in reasonable agreement with a soft X-ray image from Yohkoh, as shown in the top panel in 
Figure 11. At least the nonlinear force- free model seems to agree better with the X-ray image than 
a linear forcc-frcc and a potential field model shown in the center and bottom panel, respectively. 
Regnier et al. (2002), also using the Grad-Rubin approach, studied the non-flaring active region 
AR 8151 in February 1998 and found that the available free magnetic energy was not high enough 
to power a flare. These results are consistent which the observation in the sense that nonlinear 
force-free field lines reasonably agree with coronal observations and a consistent flaring activity: 
The particular active regions flared (not flared) when the free magnetic energy computed with 
NLFFF-codes was high enough (too low). A decreasing free magnetic energy during flares has 
been confirmed in several studies. Thalmann and Wiegelmann (2008) and Thalmann et al. (2008), 
using the optimization approach, found that the force-free energy before a small C-class flare 
(observed in active region 10960 on 2007 June 7) was 5% higher than the potential field energy. 
Before a large M-class flare (observed in active region NOAA 10540 in January 2004) the force- 
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free energy exceeded the potential field energy by 60%. In a statistic study based on 75 samples 
extrapolate with the optimization approach, Jing et al. (2010) found a positive correlation between 
free magnetic energy and the X-ray flaring index. It seems that we can trust that there is a relation 
between computed free energy and flaring activity, whereas the results of Section 6.6) indicate that 
one might not fully trust in the exact numbers of magnetic energies computed with one NLFFF-code 
only. Recently Gilchrist et al. (2012) pointed out that uncertainties in the vector magnctograms 
likely result in underestimating the computed magnetic energy. NLFFF-codes are, however, a 
strong tool to guide the investigation of coronal features. Regnier and Amari (2004), Valori et al. 
(2012), and Sun et al. (2012) applied the Grad-Rubin, MHD-relaxation and optimization approach, 
respectively and found at least qualitatively a good agreement of NLFFF-models with observed 
sigmoid or serpentine structures. 




Figure 11: Yohkoh soft X-ray image overlaid with magnetic field lines from different models: (top) 
nonlinear force-free, (center) linear force-free, and (bottom) potential fields. The original figure 
was published as Figure 8 in Bleybel et al. (2002). 
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7 Summary and Discussion 



In this review we tried to give an overview of force-free magnetic fields, particularly model as- 
sumptions, which are important for understanding the physics of the solar corona. While the 
underlying mathematical equations describe stationary states and look relatively simple, solving 
them is by no means a trivial problem because of the nonlinear nature of the problem. Exact 
solutions are only available for further simplifications, like linearizing the equations or to restrict 
to 1D/2D for the nonlinear case. For force-free configurations in 3D, we know that (for given flux 
distributions in the photosphere) the magnetic field energy is bounded from below by a potential 
field. An upper-limit for the energy is more difficult to obtain. While the Aly-Sturrock conjecture 
(Section 5.3) claims that the upper limit is for the configurations with all magnetic field lines open, 
Choe and Cheng (2002) constructed solutions with energies above the Aly-Sturrock limit. These 
configurations contain discontinuities and the debate of the validity of the Aly-Sturrock limit is 
ongoing (Hu, 2004; Wolfson et at, 2012). 

For practical computations of the 3D-ficld in the solar corona, one has to use numerical compu- 
tations and several codes have been developed, compared, and applied. As input these codes require 
measurements of the magnetic field vector in the solar photosphere. However, the transverse field 
component contains an ambiguity in the azimuth, which has to be resolved before the data can 
be used for coronal magnetic field modeling. The accuracy of photospheric measurements is lower 
for the transverse field component compared with the line-of-sight field, and in weak field regions 
measurements and azimuth ambiguity removal are less trustworthy. Consequently the majority 
of coronal force-free field models are carried out in active regions, although methods for full-disk 
computations have been developed too. A further complication of using photospheric measure- 
ments as the boundary condition for force-free extrapolations is that the photospheric magnetic 
field is not necessarily consistent with the force-free assumption. Possible solutions are to use only 
the vertical magnetic field and the vertical electric current as boundary conditions, as done for the 
Grad-Rubin approach, to preprocess the photospheric measurements with the aim to make them 
compatible with force-free and other physical requirements, or to allow changes of the transverse 
magnetic field during the iteration of a force-free field. The latter approach has been implemented 
in the optimization approach and allows us to take measurement errors into account. 

A major source for future research on force-free fields is SDO/HMI, which measures the photo- 
spheric magnetic field vector on the full disk, which in principle allows us to compute global coronal 
models as well as selecting appropriate isolated active regions with a sufficiently large field-of-view. 
Research on Stokes inversion, azimuth ambiguity removal, and force-free modeling for SDO/HMI 
data is ongoing. Another important aspect on coronal modeling is the comparison of force-free 
models as extrapolated from photospheric measurements with coronal images as observed, for ex- 
ample, with the Atmospheric Imaging Assembly (AIA; Lemen et al, 2012) on SDO. On the one 
hand, such a comparison is important to validate the models (see DeRosa et at, 2009, for details), 
and, on the other hand, the 3D models help to interpret the observations. With the 3D structure 
of magnetic loops from the models in hand, one has important tools for modeling of plasma loops, 
and may gain understanding of coronal heating and plasma flows along the loops. Further steps on 
the research of eruptive phenomena like flares and CMEs are planned with time-dependent MHD 
simulations. Force-free models are planned to be used as initial equilibria, which are disturbed 
by photospheric plasma flows (which can be deduced, e.g., from measurements with SDO/HMI). 
The temporal evolution and the potential occurrence of eruptions can be investigated with ideal 
or resistive MHD simulations in comparison with observations. Questions are if or to which extent 
the configurations remain approximately force-free during eruptions, the role of thin current sheets 
and discontinuities, and the energy and helicity content. We aim to report about the progress in 
these aspects in an update of this review in due time. 
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